// Copyright 1998-2016 Glenn McIntosh
// licensed under the GNU General Public Licence version 3
#pragma once

/** @file smooth.h
	Statistical distribution functions.
	*/

// include files
#include "linear.h"
#include <vector>

namespace math
{
using real = double;
using Vector = std::vector<real>;

/** Calculate Savitsky-Golay smoothing filter
	@param N number of points to be included (left and right)
	@param ld order of derivative (0 for smooth only)
	@param m order of smoothing polynomial
	@return calculated convolution filter
	*/
template<size_t N, int m> Array<N*2+1> savitsky_golay_filter(int ld)
{
	int np{N};

	// check arguments
	if (ld > m || np+np < m)
		throw std::invalid_argument("array too short");

	// initialize least squares matrix
	Matrix<m+1, m+1> a;
	for (int k = -np; k <= np; k++)
	{
		real exp = 1.0;
		for (int ipj = 0; ipj < m; ipj++)
		{
			for (int imj = 0; imj <= ipj; imj++)
				a[imj][ipj-imj] += exp;
			exp *= k;
		}
		for (int ipj = m; ipj <= m*2; ipj++)
		{
			for (int imj = ipj-m; imj <= m; imj++)
				a[imj][ipj-imj] += exp;
			exp *= k;
		}
	}

	// solve least squares fit
	ArrayIndex<m+1> indx;
	real d;
	luDecompose(a, indx, d);

	// solve for particular derivative
	Array<m+1> b{};
	b[ld]=1.0;
	luBackSubstitute(a, indx, b);

	// calculate coefficients
	Array<N*2+1> c;
	for (int k = 0; k <= np*2; k++)
	{
		real kk = k <= np ? k : k-(np+np+1);
		real exp = 1.0, sum = 0.0;
		for (int mm = 0; mm < m+1; mm++)
		{
			sum += b[mm] * exp;
			exp *= kk;
		}
		c[k] = sum;
	}
	return c;
}

/** Smooth a vector using an FIR filter.
	@param data data points to be smoothed in place (sequenced)
	@param respns calculated convolution filter
	@param wrap end-point algorithm (0 to extend end point, 1 to wrap around)
	*/
void fir_smooth(Vector &data, Vector &respns, bool wrap);

/** Smooth a vector using a localized Gaussian average.
	@param x independent axis (ordered and roughly equally spaced)
	@param data data points to be smoothed in place
	@param b standard deviation width of averaging weights
	*/
void gaussian_smooth(const Vector &x, Vector &data, real b); // throw(Error);
}
